# Automatic network construction and module detection without trains
library(optparse)
library(WGCNA)

option_list = list(
  make_option(c("-f", "--file"), type="character", default=NULL,
              help="exprFilePath", metavar="character"),
  make_option(c("-o", "--output"), type="character", default="out.txt",
              help="output file name [default= %default]")
  make_option(c("-o", "--output"), type="character", default="out.txt",
              help="output file name [default= %default]")
)

# exprFilePath = "/Users/zyl/Workspace/hzau/MOIN-Extension-Github/WGCNA-Extension/FemaleLiver-Data/datExpr.csv"

exprFilePath = args[1]
outputPath = args[2]

datExpr = read.csv(exprFilePath, stringsAsFactors = FALSE, row.names = 1);
nGenes = ncol(datExpr);
nSamples = nrow(datExpr);

net = blockwiseModules(datExpr, power = 6,
                       TOMType = "unsigned", minModuleSize = 30,
                       reassignThreshold = 0, mergeCutHeight = 0.25,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs = TRUE,
                       saveTOMFileBase = "/Users/zyl/Workspace/hzau/MOIN-Extension-Github/WGCNA-Extension/FemaleLiver-Data/TOM",
                       verbose = 3)

# save net data
save(net, file="network.RData")


# Visualizatioin network modules
png(
  filename = "/Users/zyl/Workspace/hzau/MOIN-Extension-Github/WGCNA-Extension/FemaleLiver-Data/network_module.png",
  width = 13,
  height = 9,
  units = "in",
  res = 300
)
mergedColors = labels2colors(net$colors)
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)

dev.off()


# Visualizatioin network heatmap
png(
  filename = "/Users/zyl/Workspace/hzau/MOIN-Extension-Github/WGCNA-Extension/FemaleLiver-Data/network_heatmap.png",
  width = 13,
  height = 9,
  units = "in",
  res = 300
)

moduleColors = labels2colors(net$colors)
dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = 6);


# # 所有基因
# plotTOM = dissTOM^7;
# diag(plotTOM) = NA;
# TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")

# 随机选择部分基因
nSelect = ceiling(nGenes / 10)
set.seed(10);
select = sample(nGenes, size = nSelect);
selectTOM = dissTOM[select, select];
selectTree = hclust(as.dist(selectTOM), method = "average")
selectColors = moduleColors[select];
plotDiss = selectTOM^7;
diag(plotDiss) = NA;
TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")


dev.off()


# Visualizing the network of eigengenes
png(
  filename = "/Users/zyl/Workspace/hzau/MOIN-Extension-Github/WGCNA-Extension/FemaleLiver-Data/network_eigengenes.png",
  width = 8,
  height = 8,
  units = "in",
  res = 300
)
MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes
MET = orderMEs(MEs)
# Plot the relationships among the eigengenes and the trait
par(cex = 0.9)
plotEigengeneNetworks(MET, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, xLabelsAngle
                      = 90)
